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We introduce a new model for an earthquake fault System that is composed of non-interacting 
simple lattice models with different levels of damage denoted by q. The undamaged lattice models 
(q — 0) have Gutenberg- Richter scaling with a cumulative exponent ß — 1/2, whereas the damaged 
models do not have well defined scaling. However, if we consider the "fault System" consisting 
of all models, damaged and undamaged, we get excellent scaling with the exponent depending on 
the relative frequency with which faults with a particular amount of damage occur in the fault 
System. This paradigm combines the idea that Gutenberg-Richter scaling is associated with an 
underlying critical point with the notion that the structure of a fault System also affects the Statistical 
distribution of earthquakes. In addition, it provides a framework in which the Variation, from one 
tectonic region to another, of the scaling exponent, or fe-value, can be understood. 

PACS numbers: 



It has been known for over half a Century that the 
frequency-magnitude distribution of seismic events oc- 
curring on a fault System obeys the Gutenbcrg-Richtcr 
(GR) scaling relation [3, 0] which states that magnitudes 
arc distributcd according to 
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where \x is the local or Richter magnitude H . Today, 
the Richter magnitude has been superseded by the more 
physically significant moment magnitude, ß w , [4| which is 
related to the average slip, Am, over the rupture surfacc 
by 
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A is the area of and G is the shear modulus of the rupture 
surface, and M is the so-called seismic moment. The 
exponential distribution in the magnitude then becomes 
a power-law distribution in the moment, transforming 
Eq. [TJ into 
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with ß = 
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The distribution Nm is the total number of events with 
moment greater than or equal to M [l], Q . The inverse 
characteristic magnitude, b, in Eq. Q]is distributcd in the 
narrow ränge 0.75 <6< 1.2 This ränge of b values 

translates into a ränge 0.5 < ß < 0.8 for the exponent of 
the moment distribution. 

The GR distribution in Eqs. [1] and E] has been the sub- 
ject of considerable investigation by seismologists as well 
as physicists and geo-physicists who have associated the 



scaling with an underlying critical point. This apparent 
critical phenomena has led to the creation of several sim- 
ple models that exhibit a GR-like distribution. The first, 
proposed by Burridge and Knopoff Q , consists of blocks 
connected by linear Springs resting on a rough surfacc 
and is evolved according to Newton's Laws. Using the 
Mohr-Coulomb friction law and ignoring inertial effects, 
Rundle and Jackson 0] introduced a cellular automaton 
version of the Burridge-Knopoff model which was rein- 
troduced in lattice form by Olami et al. (OFC) @. In 
addition to the association with critical points, Bäk et 
al. [IÖ| and Feder and Feder [ll| proposed that GR scal- 



ing is a form of seif organized criticality. Klein et al. [12 1 
proposed that the scaling is associated with fluetuations 
about a spinodal critical point. All of these models share 
at least one shortcoming: they are best understood as de- 
scribing individual faults rather than fault Systems while 
it is the latter which are observed to scale according to 
Eqs. Q] and [3] Seismic events on Single faults may or may 
not be described by the GR distribution, depending on 
the fault, as we will show later. 

In this letter, we propose a different approach that 
models a fault System as a collcction of non-interacting 
OFC models or "faults" with varying degrees of dam- 
age. The frequency of events on individual faults does 
not, in gcneral, scale as a power-law in the event's mo- 
ment, however, when a collection of model faults with 
varying degrees of damage are assembled into a "fault 
System", the aggregate distribution of moments is well 
described by the GR scaling relation. The OFC model 
consists of a regulär lattice of linear size L with open 
boundary conditions and a stress failure threshold, a F , 
and a residual stress, a R , at each vertex (site). Initially, 
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the sites are assigned a random stress, a R < aj < a F , 
whcrc thc index j labcls the position of the site on the 
lattice. The model is evolved via a plate update in 
which stress is simultaneously added to each site until 
the stress on the site nearest its failure threshold reaches 
a F . This loading method is the so-called zero velocity 
limit Q. When a site fails its stress is reduced to a R 
and (1 — — cr R )/ar 2 is distributed to each of its ar 2 
neighbors where a is a dimensionless number depending 
on the shape of the interaction, r is the ränge of the in- 
teraction, and < a < 1 is the dissipation parameter. 
Failure of a site may cause ncighboring sites to fail, re- 
sulting in an carthquakc or avalanche. We continue to 
scarch the lattice until crj < <j f V j, at which point the 
lattice is reloadcd as dcscribed above. 

We add a small amount of noise to the model by re- 
sctting each site to a R ± rj, where 77 is drawn from a 
flat distribution with mean a R . For this paper, we take 
a F = 2, a R = 1 ±0.025, L = 512, and stress from a failed 
site is redistributed within a circle of radius r = 20. We 
choose r = 20 to mimic the long-range stress transfer 
in fault Systems 



We simulate the System for 
4 x 10 6 plate Updates to remove any transient effccts be- 
töre recording the size (number of failed sites) of each of 
the next 10 7 events. 

For nearest neighbor stress transfer (r — 1), there is no 
apparent GR scaling [l5 - 17 1 , however if we take r>l, 
the System becomes near meanfield 15[ and there ap- 
pears to be GR or Cluster scaling associated with a spin- 
odal critical point 12, Hj]. In this limit, the number of 
events, n s , of size s£ [s , s + ds] scales as 
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with r = 3/2. Here, no is a measure of the seismic ac- 
tivity of the fault and Ah is a measure of the distance 
from the spinodal [15l . Il8l ] . To make contact with obser- 
vations on real faults, we note that for r>l, the size, s 
in Eq. [4l is proportional to the moment (l5| . To obtain 
the cumulative distribution, Nm, we simply integrate 
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We now modify the model by randomly eliminating a 
fraction of sites, q [ljjj]- These eliminated or dead sites 
do not hold stress, hence, when stress is transferred to a 
dead site from a failing neighbor, it is dissipated. In real 
faults, this can be caused by different levels of fracture or 
gouge, which are known to vary over a fault System [2Ö| . 
We do this for many lattices with < q < 1. As can 
be seen from Fig[TJ the larger the value of q the less the 
distribution resembles a simple power-law. This is con- 
sistent with the idea that the OFC model in the zero 
velocity limit is only critical with a zero dissipation con- 
stant [21[. Nevertheless, the data is still well describcd 
by Eq. 31 but with q-dependent Ah and ng. 
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FIG. 1: (color online) The frequency-size distribution for a 
model fault for various values of q. For q = (i.e. no dam- 
age) we observe scaling up to the linear System size. As q is 
increased, the distribution is no longer described by a power- 
law, rather, large events are exponentially suppressed aecord- 
ing to Eq. [4] The solid black line has slope —3/2, the expo- 
nent predicted for the scaling on undamaged faults. To see 
the distributions clearly, we do not plot n s for every value of 
s. Instead, we plot the distribution at evenly spaced values 
of s in log-space. All distributions contain the same number 
of events. 



To determinc the dependence of Ah and no on 5, we 
consider the theoretical description of the model as de- 
rived by Klein et al. 



12L Il5l ] from which one can argue 



for thc dependence 
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where r = 3/2, as above, and is a q-independent con- 
stant related to the average seismic activity on the fault. 
We check this g-dependence in two ways. First, for var- 
ious valucs of q, wc fit n s , using a least-squarcs, non- 
linear, weighted fit, for Ah and no to Eq.[6j Additionally, 
we introduce a scaling variable z = q 2 s such that plots 
of n z (1 — q)/q 3 should collapsc to a Single curve for all 
values of q. The data collapse in Fig. [2] suggests we have 
the correct scaling form for n s . 

If this mechanism is rcsponsiblc for thc frequency- 
magnitude distribution on single faults, then we expect 
the data obtained from real faults to vary from fault to 
fault but, this Variation should be captured by a Single 
parameter in an otherwise universal distribution, as in 
Eq. 01 Using scaling arguments, one can rewrite Eq. [Hin 
terms of the slip area, A, introduced in Eq. [2 as 
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FIG. 2: (color online) The frequency-size distribution plotted 
in terms of the scaling variables. We explicitly account for the 
g-dependence in n s by plotting (1 — q)/q z n z versus z = q 2 s. 
Having removed the g-dependence in these distributions, the 
data collapses to a Single curve for all values of q. To see 
the collapse clearly, we plot the distribution at evenly spaced 
values of s in log-space. All distributions contain the same 
number of events. 



In Fig. [3j we plot the cumulativc number of events with 
area A on three faults from the southern California Sys- 
tem. As can be seen from the figure, the data is consistent 
with Eq. [3 

To obtain the frequency-size distribution for a fault 
system, we must integrate Eq. [6]over all values of q, that 
is, 
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where the weighting factor, D q , is the fraction of faults 
with damage q 6 [q , q + dq] . For equal weighting 
(D q = 1) we find n s ~ l/s 2 which implies a cumula- 
tive exponent ß = f — 1 = 1 (c.f. Eq. [5]). In Fig. |H 
we plot (purple crosses) the number of events of a given 
size versus s. As noted above, ß is observed to vary be- 
tween 0.5 and 0.8 in real fault Systems 0, [öjj, however, 
in real fault Systems we do not know D q . If we associate 
the damage with micro-cracks, which are know to have 
a fractal spatial distribution [20J, it is reasonablc to as- 
sume a scale-free or power-law distribution for D q ~ q~ x , 
where the exponent varies depending on the fault System. 
Thercfore, the frequency-size distribution in Eq.[8]can be 
written in terms of the incomplcte gamma funetion. If 
we perform an asymptotic expansion in the limit s>l, 
we find 
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FIG. 3: (color online) The number of events with magnitude 
greater than or equal to the minimum ma gnit ude of complete- 
ness [13] versus the area of each event [23l | occurring from 
1980 - 2008 within a 20 km swatch on either side of the San 
Jacinto fault (green crosses), Ft. Tejon segment of the San 
Andreas Fault (blue diamonds), and the Creeping section of 
the San Andreas Fault (red triangles). To see the data clearly, 
the various faults have been offset by factors of ten. The fits 
(solid colored lines) are to Eq. [7] and the solid black line has 
unit slope. 



where T(z) is the gamma funetion. To leading order, we 



find 



with f = 2 — x/2, and so if we wish to 



obtain a 6-value of unity, we simply solve 
,3(2-0;) 
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with b = 1, so x = 2/3. Furthermorc, it is clear from 
Eq.[TÜ]that a ränge of b values can be obtained by varying 
the exponent in D q [24[ as shown in Fig [4] In Fig [5l 
we plot the the scaling exponents f and b versus the 
exponent, x, of the damage distribution, D q . The solid 
red line is the rclationship predicted by Eq. [9] while the 
black exes represents the value of f determhied by a linear 
fit of the data to the form N/s T for N and f. To get 
better agreement between this fit value and our predicted 
value, we perform a two parameter non-lincar fit to the 
form N/s f {l+~ f /s 1 / 2 ) where 7 = r[(2-a;)/2]/r[(l-a;)/2] 
is not a fit parameter. The value of f as determined from 
this fit is shown in blue crosses. 

In summary, we have introduced a new paradigm for 
understanding GR scaling which links scaling associated 
with a (spinodal) critical point on a singlc fault with GR 
scaling on a fault System. The GR scaling on a fault 
System is affected by both the degree of damage on each 
fault (q), and the distribution of damage over the fault 
System (D q ). This paradigm offers an explanation for the 
Observation of GR scaling on fault Systems comprised of 
individual faults: none of which produce the same event 
magnitude statistics as the aggregate System. It also pro- 
vides a framework in which one can understand the ob- 
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served Variation of thc GR exponcnt, or 6- value, from one 
rcgion to anothcr. In addition, since the various faults 
that comprise the fault System are non-interacting by 
construction, our model provides a mechanism in which 
scaling can arise on fault System length scales that are 
ten to one-hundred times larger than any length scale 
associated with the critical phenomena occurring on a 
Single fault. 
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FIG. 4: (color online) The cumulative frequency-size distribu- 
tion for a model fault System for various damage distributions, 
D q ~ q~ x . By varying the exponent x, the GR distribution 
scales with different exponents f. The solid black lines have 
slope —19/12 and — 2, the exponents predicted by Eq. |Ü]for 
x — 5/6 and x = 0, respectively. To see the data clearly, 
we plot the distribution at evenly spaced values of s is log- 
space and we offset different values of x by scaling the data 
by factors of 10. All distributions contain the same number 
of cvcnts. 
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FIG. 5: (color online) The scaling exponents, f and b (c.f 
Eqs. Q] and [H] respectively), versus the damage distribution 
exponcnt, x. The black exes represent the value of f as de- 
termined by Atting the data to a Single power-law. The blue 
crosses are the value of f detcrmined by Atting the data to a 
power-law plus the leading asymptotic correction to the scal- 
ing form. The solid red line is the value of f predicted by 
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